#pragma once

#include "Tensor.h"
#include "Polynomial.h"

/**
 * @file Integrator.h
 * @author xingyifan
 * @date 2022-08-31 22:34
 *
 * @description: 网格容器
 */

template <int Dim>
using Grid = Tensor<Dim>;

template <int Dim>
using GridX = std::vector<Grid<Dim>>;

/**
 * @file Integrator.h
 * @author xingyifan
 * @date 2022-08-31 22:35
 *
 * @description: 用于离散化的区域
 */

template <int Dim>
class Domain
{
public:
    // 默认生成单元区域
    Domain() : width_(1)
    {
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = 0;
        }
    }
    Domain(Real *l, Real width) : width_(width)
    {
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = l[i];
        }
    }
    Domain(std::initializer_list<Real> l, Real width) : width_(width)
    {
        auto it = l.begin();
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = *it++;
        }
    }
    Domain(const Domain &other) : width_(other.width_)
    {
        for (int i = 0; i < Dim; i++)
        {
            lower_[i] = other.lower_[i];
        }
    }

    Real width() const { return width_; }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:35
     *
     * @description: 在 d 维上移动区域
     */

    void move(int d, int dist) { lower_[d] += dist * width_; }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:36
     *
     * @description: 计算区域体积
     */

    Real volume() const { return pow(width_, Dim); }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:36
     *
     * @description: t 当前指标对应的格点通过 Domain 转换为区域中对应的格点
     */

    Vector gridPoint(const Grid<Dim> &t) const
    {
        Real h = width_ / (t.base() - 1);

        // 获得角点
        Vector point(Dim);
        for (int i = 0; i < Dim; ++i)
        {
            point[i] = lower_[i] + t.index(i) * h;
        }
        return point;
    }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:36
     *
     * @description: 用于有限体积法，t 当前指标对应的格点通过 Domain 转换为区域中对应的单元角点
     */

    Vector cell(const Grid<Dim> &t) const
    {
        Real h = width_ / t.base();

        // 获得角点
        Vector point(Dim);
        for (int i = 0; i < Dim; ++i)
        {
            point[i] = lower_[i] + t.index(i) * h;
        }
        return point;
    }

    friend std::ostream &operator<<(std::ostream &out, const Domain &d)
    {
        for (int i = 0; i < Dim; i++)
        {
            out << "[" << d.lower_[i] << "," << d.lower_[i] + d.width_ << "]" << std::endl;
        }
        out << std::endl;
        return out;
    }

private:
    const Real width_; // 宽度
    Real lower_[Dim];  // 下界
};

/**
 * @file Integrator.h
 * @author xingyifan
 * @date 2022-08-31 22:36
 *
 * @description: 多元数值积分 D 阶，Dim 维
 */

template <int D, int Dim>
class Integrator
{
public:
    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:37
     *
     * @description: 在构造函数中初始化权重
     */

    Integrator()
    {
        // 计算权的分布
        Tensor<1> baseWeight(D + 1);

        // 使用插值多项式积分计算权
        std::vector<Real> x(D + 1);
        for (int i = 0; i < D + 1; i++)
        {
            x[i] = i;
        }
        for (int i = 0; i < D + 1; i++)
        {
            std::vector<Real> f(D + 1, 0);
            f[i] = 1;
            baseWeight[i] = Polynomial::interpolate(f, x).intergral(0, D);
        }

        // 根据维数取权张量的积
        weight_ = TensorOp::power<Dim>(baseWeight);
    }

public:
    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:37
     *
     * @description: 向量函数积分：FuncX 接收 Dim 元向量，输出 n 元向量
     */

    Vector NewtonCotes(Domain<Dim> &domain, FuncX f)
    {
        // 分割宽度
        Real h = domain.width() / D;

        // 离散 Dim 维区域为网格
        Grid<Dim> grid(D + 1);

        // 获得输出向量的长度
        int n = f(Vector(Dim, 0)).size();
        Vector I(n);

        // 逐分量积分
        for (int i = 0; i < n; i++)
        {
            // 填充区域
            do
            {
                Vector y = f(domain.gridPoint(grid));
                // 获得对应区域中的点计算函数值
                grid.value() = y[i];
            } while (grid++);

            // 获得单元值并加权平均
            Tensor<Dim> value = grid.cell();
            I[i] = TensorOp::dot(weight_, grid) * pow(h, Dim);
        }

        return I;
    }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:37
     *
     * @description: 计算对应指标区域上的体元平均，dist 用于移动单元，d 表示在 d 维上移动
     */

    Vector CellAverage(Domain<Dim> &domain, Grid<Dim> &index, FuncX f, int d = 0, int dist = 0)
    {
        // 获得指标对应的单元角点
        Vector point = domain.cell(index);

        // 单元长度
        Real h = domain.width() / index.width();

        // 用这个点作为角点，构造积分区域
        Domain<Dim> cell(point.smartptr().get(), h);

        // 移动单元
        cell.move(d, dist);

        return NewtonCotes(cell, f) / cell.volume();
    }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:38
     *
     * @description: 计算对应指标位置在 d 维的左/右面元的平均
     */

    Vector FaceAverage(Domain<Dim> &domain, Grid<Dim> &index, int d, int dist, FuncX f)
    {
        // 获得指标对应的单元角点
        Vector point = domain.cell(index);

        // 单元长度
        Real h = domain.width() / index.width();

        // 构造积分区域
        Domain<Dim> cell(point.smartptr().get(), h);

        // 根据方向移动区域，dist 中 -1,1 表示左右面，如果增加 -2,2，则表示移动一个单元后的左右面
        if (dist > 0)
        {
            cell.move(d, dist - 2);
        }
        else
        {
            cell.move(d, dist - 1);
        }

        // 存放积分结果
        Vector res = NewtonCotes(cell, f) / cell.volume() * (-1.0 / 12);

        // 回到原处 phi
        cell.move(d, 1);
        res = res + NewtonCotes(cell, f) / cell.volume() * (7.0 / 12);

        // phi_ed
        cell.move(d, 1);
        res = res + NewtonCotes(cell, f) / cell.volume() * (7.0 / 12);

        // phi_2ed
        cell.move(d, 1);
        res = res + NewtonCotes(cell, f) / cell.volume() * (-1.0 / 12);

        return res;
    }

    /**
     * @file Integrator.h
     * @author xingyifan
     * @date 2022-08-31 22:39
     *
     * @description: 用于体积分填充
     */

    void fill(Domain<Dim> &domain, Grid<Dim> &index, FuncX f)
    {
        // 遍历每一格点填充
        do
        {
            index.value() = (Real)CellAverage(domain, index, f);
        } while (index++);
    }

    void fill(Domain<Dim> &domain, std::vector<Grid<Dim>> &index, FuncX f)
    {
        // 输出向量的规模
        int n = index.size();

        // 检测退出
        bool check = true;

        // 遍历每一格点填充
        do
        {
            Vector F = CellAverage(domain, index[0], f);
            for (int i = 0; i < n; i++)
            {
                index[i].value() = F[i];
                check = index[i]++;
            }
        } while (check);
    }

public:
    Tensor<Dim> weight_; // 积分权
};